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Abstract 

In a recent Letter [Phys. Rev. Lett. 103, 097403 (2009)], we outlined a computational method 
to calculate the optical properties of structures with a spatially nonlocal dielectric function. In this 
Article, we detail the full method, and verify it against analytical results for cylindrical nanowires. 
Then, as examples of our method, we calculate the optical properties of Au nanostructures in one, 
two, and three dimensions. We first calculate the transmission, reflection, and absorption spec- 
tra of thin films. Because of their simplicity, these systems demonstrate clearly the longitudinal 
(or volume) plasmons characteristic of nonlocal effects, which result in anomalous absorption and 
plasmon blueshifting. We then study the optical properties of spherical nanoparticles, which also 
exhibit such nonlocal effects. Finally, we compare the maximum and average electric field enhance- 
ments around nanowires of various shapes to local theory predictions. We demonstrate that when 
nonlocal effects are included, significant decreases in such properties can occur. 

PACS numbers: 78.67.Lt, 78.67.Uh, 78.20.-e, 77.22.Ch 
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I. INTRODUCTION 



Interest in the optical properties of metallic nanostructures has been steadily increas- 
ing as experimental techniques for their fabrication and investigation have become more 
sophisticated^. One of the main driving forces of this is their potential utility in sens- 
ing, photonic, and optoelectronics applications^^. However, there can also be interesting 
fundamental issues to consider, particularly as very small length scales are approached (ap- 
proximately less than 10 nm). In this limit, quantum mechanical effects can lead to unusual 
optical properties relative to predictions based on classical electrodynamics applied with 
bulk, local dielectric values for the metaP. In isolated spherical nanoparticles, for example, 
localized surface plasmon resonances (LSPRs) are found to be blueshifted relative to Mie 
theory predictions^, and in thin metal films, anomalous absorption is observed^. 

Roughly speaking, when light interacts with a structure of size d (e.g., a nanoparticle size 
or junction gap distance), wavevector components k, which are related to the momentum 
p by p = hk, where % is Planck's constant, are generated with magnitude k = 2it/d. This, 
in turn, imparts an energy of E — (hk) 2 /2m e , where m e is the mass of an electron, to 
(relatively) free electrons in the metal. For small d, these energies can correspond to the 
optical range (1-6 eV). This analysis suggests that such effects should come into play for 
d less than approximately 2 nm. In metals, however, somewhat larger d values also exhibit 
these effects, because electrons in motion at the Fermi velocity can be excited by the same 
energy with a smaller momentum increase, due to dispersion effects. 

A full quantum mechanical treatment of such structures would of course be best, but this 
is not practical for these sizes. However, it is possible to incorporate some "quantum effects" 
within classical electrodynamics via use of a different dielectric model than that for the bulk 
metal. At least four such effects can be addressed in this way: electron scattering, electron 
spill-out, quantum size-effects, and spatial nonlocality of the material polarization. The 
additional losses due to increased electron scattering at the metal surface can be described 
by a size-dependent damping ternP, which will effectively broaden spectral peaks^, as we 
consider below. The electron spill-out from the metal into the medium, due to the electron 
density varying smoothly, can partially be accounted for by a dielectric layer model. The 
effect of this is varied, and depends on a number of details, including the surface chemistry 
of the structure^ and its local dielectric environment. Quantum size effects due to discrete 
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electronic energy levels can lead to a size and shape-dependent conductivity. At least for 
metal filmd-^, this quantity is reduced relative to the bulk, and exhibits peaks for certain 
film thicknesses. Such effects can be incorporated directly into classical calculations for 
some simple systems, based on rigorous theory. Although, because this effect and electron 
spill-out are both highly dependent on system specifics, incorporating them into a general 
framework is not straightforward. Therefore, they will not be considered in this work. The 
fourth quantum effect, and the one that is of main interest herein, is the need of a dielectric 
model which considers that the material polarization at one point in space depends not only 
on the local electric field, but also that in its neighborhood^. 

In classical electrodynamics, materials are described through a dielectric function e that 
relates the electric displacement field D (proportional to both the incident field and material 
polarization) to the electric field E at a given frequency of light u. This relationship is usually 
assumed to be local in space (i.e., the polarizability of the material at a point x only depends 
on E at x). However, in the more general case 

D(x, to) = Eq J c/x'e(x, x', w)E(x', to) , (1) 

where e(x, x', u) is a spatially-dependent (nonlocal) and frequency-dispersive relative dielec- 
tric function. In a homogeneous environment (an approximation which we make for the 
finite, arbitrarily shaped structures considered herein), e(x, x',o;) only spatially depends on 
|x — x'|. Therefore, in k-space, Eq. ([T]) is more simply expressed as 

D(k,w) = e e{k,u)E{k,u) . (2) 

Since the first formulation of nonlocal electromagnetics^, applications of k-dependent 
dielectric functions have remained limited to simple systems, such as spherical structure d 13 1 14 1 
or aggregates thereoP^^, and planar surfaces^. Nonetheless, this k-dependence has been 
found experimentally^ and proven theoretically * 13 * 19 * to have important consequences. For 
example, such dependence is responsible for the aforementioned anomalous absorption and 
LSPR blueshifting. 

In a recent Letter^, we outlined a powerful, yet simple method by which the optical 
properties of arbitrarily shaped structures with a nonlocal dielectric function can easily 
be calculated. This was done by deriving an equation of motion for the current associ- 
ated with the hydrodynamic Drude modeP^, which we solved within the framework of the 
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finite-difference time-domain (FDTD) method 22 . The advantage of this approach is that it 
can describe the dynamical optical response of structures that are too large to treat using 
quantum mechanics, yet small enough such that the application of local continuum electro- 
dynamics becomes questionable. In this Article, we expand on that work and detail the full 
method. We first verify our results against analytical ones for the cylindrical Au nanowires 
that we considered in our previous Letter^. Then, as new examples, we calculate the optical 
properties of one, two, and three dimensional Au nanostructures. 



II. THEORETICAL APPROACH 
A. Formulation of the Method 



The interaction of light with matter in the classical continuum limit (i.e., many hundreds 
of atoms or more) is described by Maxwell's equations, 

^D(x,t)+J(x,t) = VxH(x,t) (3) 

^B(x,t) = -VxE(x,t) (4) 

V-D(x,t)=p (5) 

V-B(x,t) = (6) 

where H(x, t) and B(x, t) are the auxiliary magnetic field and the magnetic field, respec- 
tively, while J(x, t) and p are external current and charge densities. Except for the most 
simple systems, such as spheres or metal films, analytical solutions or simplifying approx- 
imations to Eqs. (J3]) - Q do not exist. Therefore, computational methods are often used 
to solve them, one of the most popular being FDTE^l. For dynamical fields, Eqs. ^ and 
Q are explicitly solved, while Eqs. ^ and (|6| are considered initial conditions that should 
remain satisfied for all time. 

However, before Maxwell's equations can be solved, an explicit form for e(k, u) in the 
constitutive relationship between D(k, lo) and E(k, uj), Eq. ([2]), must be specified. Note 
that Eqs. @ - Q are in terms of x and t, but material properties are often dependent 
on k and u, which are related to the former via Fourier transform. Also note that we 
assume that there are no magnetic materials present, and thus the magnetic field constitutive 
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relationship is B(x, oj) = /xoH(x, u), where //q is the vacuum permeability. Returning to the 
current discussion, the dielectric function of a metal like Au is well-described in the classical 
continuum limit by three separate components, 

e(k, U) = Boo + SinterM + S in tra(k, u) , (7) 

the value as u — > oo, e^, a contribution from <i-band to sp-band (conduction band) inter- 
band electron transitions, Sinter (w), and a contribution due to sp-band electron excitations, 
£intra(k, of). The notation in Eq. (|7| highlights the k and u dependencies. 

Winter (^) can be physically described using a multipole Lorentz oscillator modeP^, 

Sinter (w) = > -3 7 , - ox x 8 

where j is an index labeling the individual d-band to sp-band electron transitions occurring 
at WLj, Asl? is the shift in relative permittivity at the transition, and 8i,j is the electron 
dephasing rate. Because there are two interband transitions in Au at optical frequencies 
(near 3 and 4 eV^), we take j = 2 in this work. 

Smtra(k, oj) is responsible for both the plasmonic optical response of metals and nonlocal 
effects. Both of these can be described by the hydrodynamic Drude modepH, which reduces 
to the local Drude expression for electron motion if k — > Cp^l, 

oj 2 

S int ra(k,U;) = ° ~o2]~2 ( 9 ) 

0J{0J + lj) — p z k 

where ojd is the plasma frequency, 7 is the collision frequency, and for a free electron gas (i.e., 
one with only kinetic energy) (3 2 = Cv^/D, where Vp is the Fermi velocity (1.39 • 10 6 m/s for 
Au), D is the dimension of the system, and C = 1 at low frequencies and 3D/(D + 2) at high 
frequencie s 2 ^ 2 1 . We note, in passing, that other analytical forms for e- mtvai (\t, oj) could also 
be used with the following approach, such as those inferred from representative quantum 
mechanical electronic structure calculation^. 

Inserting Eqs. ^ and ^ [using Eqs. ^ and ([9])] into the Maxwell- Ampere law in k-space 
for a time-harmonic field, — zwD(k, 00) = ik x H(k, oj), leads to 

- iue e oo E(k, oj) + J L j(k, oj) + J H p(k, oj) = ik x H(k, oj) , (10) 



where the Jl? (k, oj) are polarization currents associated with Eq. ([8 



and JnD(k, u>) is a nonlocal polarization current associated with Eq. ([9]), 



J H D(k,o;) 



LOT 



rE(k,w) 



(12) 



u{u + ij) - P 2 k 2 ' 

Equations of motion (i.e., partial differential equations in terms of x and t) for the currents 



in Eqs. (11) and (12) can be obtained by multiplying through each by the appropriate 
denominator and inverse Fourier transforming (zk — > V and — iu — » d/dt), 

— J L ,(x, t) + 25 Lj — J Lj (x, t) + ^J Lj (x, t) = eoAe^ulj— E(x, t) (13) 



at 
o 2 



'dt 



, !i2 J H D(x,t)+7^JHD(x,t)-/3 2 V 2 J H D(x,t) = e u 2 D -E(x,t) . (14) 



Because of the spatial derivatives in Eq. (14), the equation of motion for the hydrodynamic 
Drude model is second-order, unlike the normal Drude model that is first-order^. 

Equations (13) and (14) can be solved self-consistently with Eq. Q and the inverse 
Fourier-transformed form of Eq. (ITo| , 







e £oo^E(x,t) +J2 j LjW + Jhd(x,£) = V x H(x,t) 



(15) 



along with the requirement that Eqs. ^ and ^ are, and remain satisfied. Our implemen- 
tation of Eqs. Q and Q - (|15( using standard finite-difference techniques is outlined in 
the Appendix. 



B. Au Dielectric Function 



To model Au nanostructures, and use the approach outlined in Subsection II A , Eq. ^ 



must first be fit to the experimentally determined dielectric data of bulk Aii^. This is done 
in the limit of k — > in Eq. Q, which is valid for large structures, such as those used to 
obtain the experimental data. To make sure that the separate terms in Eq. ^ accurately 
capture the physics of the problem, it is necessary to fit Eqs. ^ and ^ over the appropriate 
energy ranges separately. Using simulated annealing, we first fit Eq. ^ (also incorporating 
£oo) over the range 1.0 - 1.8 eV, where e(0,u) is dominated by sp-band electron motion. 
Then, keeping the parameters in Eq. (|9| constant (but not £oo), the entire dielectric function 
in Eq. ^ was fit over the full range of interest, 1.0 - 6.0 eV. The resulting parameters were: 
Eoo = 3.559, co D = 8.812 eV, 7 = 0.0752 eV, A^li = 2.912, co L1 = 4.693 eV, S L1 = 1.541 eV, 
Ae L2 = 1.272, u L2 = 3.112 eV, and 5 L2 = 0.525 eV. 



A plot of calculated dielectric values against those experimentally determined is shown 
in Fig. [I] The fit is reasonably good given the simple form of Eq. ([7]). For example, features 
of the two interband transitions are captured near 3.15 and 4.30 eV, which is evident in 
lmag[e(0, co)}. Note that Co>li and cjl2 are also close to these values, as expected based on the 
discussion in Subsection II A Although, the fit is not as good as could be achieved with a 
more flexible function, such as an unrestricted fit. However, the present fitting scheme leads 
to parameters that are more physically realistic, and this is essential given that we are going 
to use these local (k = 0) parameters in the nonlocal (k ^ 0) expression. One consequence 
of this fit, for example, is that the minimum value of lmag[£:(0, u)] near 1.85 eV is not as 
small as the experimental one, which will end up giving broader plasmon near this energy 
resonances than expected. Fortunately, such differences will not play a significant role in 
the results that we present. 

It is interesting to look at the evaluation of Eq. (m) as a function of both k and u; 
Fig. [2] Note that a slice through k = gives the local dielectric data in Fig. [T] When 
(3k <C u, e(k, u) is relatively constant for a given u - i.e., it remains close to the local value. 
However, as (3k approaches oj from below, e(k, oj) quickly becomes very negative and changes 
sign rapidly as it passes through (3k w u, after which e(k, oj) is thus no longer plasmonic. 
Absorption of light by a system is related to the value of e(k, oj) and the structure under 
consideration. For example, for a small spherical particle in air, the maximum absorption 
occurs when Real[e(k, oj)] = — 2 2 - : . Figure [2] therefore indicates that in addition to the local 
absorption, additional (anomalous) absorption will occur when (3k w oj {when the rapid 
variation in Real[e(k, oj)] occurs}. 

Nonlocal effects are very prominent for small structures^, as we will demonstrate below. 
In such systems, it is necessary to consider the reduced mean free path of the sp-band 
electrons due to electron-interface scattering. As was briefly discussed in Section [I], this 
can be taken into account by using a modified collision frequency in Eq. ([9]) 8 : 7' = 7 + 
AvF/L e s, where the effective mean free electron path is L c g = 4V/S in 3D and nS/P 
in 2D, where V is the volume of the structure with surface area S with perimeter P, 
and A can be considered the proportion of electron-interface collisions that are totally 
inelastic. Such scattering can also be considered a nonlocal effectP^^. In a formal sense, A 
is related to the translational invariance at the surface, the full description of which depends 
on the geometry and morphology of the structure, its local dielectric environmentP^, and 
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FIG. 1: (color online) Fitted dielectric data for bulk Au (solid red lines), compared to that exper- 
imentally determined (open blue circles^PH. 

the dielectric function of the material, which is ultimately nonlocal in character. Although, 
the general magnitude of A can be arrived at in the local limit, and in a variety of wayiP. 
Because of these complex details, correctly choosing the value of A can be challenging, and 
large values can have a significant effectP. Therefore, for simplicity, and consistency with 
our previous Letter^, we take A = 0.1 for the calculations herein. 
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FIG. 2: (color online) |Real[e(k, w)]| of Au as a function of both k and w. Below /3/s ~ uj, e(k, w) < 0; 
and above, e(k, cj) > 0. Note that no specific value is attached to (3 2 . The condition /3k = ui is 
shown using a dashed white line. 

III. COMPUTATIONAL CONSIDERATIONS 

Computational domains were discretized using a Yee spatial lattice^, as outlined in the 
Appendix. The edges of the domains were truncated using convolutional perfectly matched 
layers^. For the calculations in Section IV, variable grid spacings were used for this dis- 



cretization, as will be discussed, as well as the low-frequency 2D value of (3 2 (for consistency 
with our previous Letter^). For the calculations on metal films in Subsection VA, grid 



9 



spacings of 0.1 nm were used for the 2-nm film, and 0.2 nm for the others, as well as the 
high-frequency 2D value of /3 2 . For the nanoparticles in Subsection VB grid spacings of 0.2 



nm were used for the 4 and 7-nm nanoparticles, and 0.5 nm for the 15-nm one. For all of 



the nanowires in Subsection VC grid spacings of 0.25 nm were used. In both of these latter 
cases, the high-frequency values of (3 2 were used. 

Optical responses were determined by calculating extinction cross sections^ (the amount 
of power absorbed or scattered relative to the incident light), by integrating the normal 
component of the Poynting vector around surfaces enclosing the particles^El. In order to 
obtain accurate Fourier-transformed fields necessary for such calculations, as well as field 
intensity profiles, incident Gaussian damped sinusoidal pulses with frequency content over 
the range of interest (1-6 eV) were introduced into the computational domains using the 
total-field-scattered-field technique^. Furthermore, all simulations were carried out to at 
least 100 fs. 

Before leaving this subsection, we mention that instabilities have been encountered in 
some 3D calculations. For example, simulations of 1.0-nm Au nanoparticles become unstable 
when grid spacings of 0.05 nm are used. However, in 2D, such instabilities do not seem to 
exist. We are currently investigating this issue. 

IV. NUMERICAL VERIFICATION 



One way to determine the accuracy of the method presented in Subsection II A and the 
Appendix is to compare computed results to obtainable analytical ones. Such comparisons 
are possible for metal films^, cylindrical wires^, and spherical particles^. In this section, 
such a comparison is made for a r = 2 nm radius cylindrical nanowire, an example that we 
considered in our previous Letter^. 

The optical responses calculated using uniform grid spacings A in both x and y of 0.2, 
0.1, or 0.05 nm are compared to the analytical result}^ in Fig. |3j A number of peaks and 
valleys are seen in all results, which correspond to the dipolar LSPR near 2.55 eV and 
anomalous absorption near 1.61, 2.75, and 3.78 eV. A discussion of these effects was given in 
Ref. EDI and will be further elaborated on below. The calculated and analytical results are 
found to agree remarkably well, both qualitatively and quantitatively, providing numerical 
verification of our method. Additionally, decreasing the grid spacing leads to significantly 
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FIG. 3: (color online) Convergence of calculated extinction cross sections to the analytical resuhP^ 
for a 2 nm radius cylindrical Au nanowire, with respect to the grid spacing A. 

better results, especially for the higher energy peaks. For example, the additional peak near 
3.78 eV converges from 3.53 to 3.68 to 3.72 eV as A is reduced from 0.2 to 0.1 to 0.05 nm, 
respectively. Such convergence is understandable, because for a given grid spacing, there is 
an uncertainty in r of ±A. Since the results appear to redshift with increasing A, we can 
infer that the nanowire radius is approximately r + A. The convergence of these effects is 
much tougher than those in local electrodynamics, where even A = 0.2 nm is sufficient (not 
shown). These results demonstrate the exquisite sensitivity of nonlocal effects to even minor 
geometric features. 

Grid spacings of 0.05 or 0.1 nm are impractical for most calculations, because of the 
resulting computational effort required. At first this appears troublesome, given that the 
nonlocal effects are so sensitive to this parameter. However, the only real downside is that 
an uncertainty in the geometry of ±A must be accepted (which is also the case in local 
electrodynamics, but is less important). This is because it is found that for a given grid 
spacing, the calculated results always fall between the analytical ones that incorporate the 
±A tolerance; see Fig. |4| For example, in the case of a r = 2 nm cylindrical nanowire, the 
calculated results with A = 0.2 are constrained by the analytical ones with r = 1.8 and 2.2 
nm. Results for A = 0.1 and 0.05 nm are shown in Fig. [4] as well, and are also consistent 
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FIG. 4: (color online) Extinction cross sections of 2 nm radius cylindrical Au nanowires, calculated 
using grid spacings of A = (top) 0.2, (middle) 0.1, and (bottom) 0.05 nm. The results show that 
the calculations (solid black lines), denoted as FDTD, are always constrained by the analytical 
results (broken lines) when tolerances for the grid spacings are considered. 



with this analysis. 
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V. APPLICATIONS 



A. Metal Films (ID Systems) 

In this subsection, we study the transmission, reflection, and absorption of thin Au films 
illuminated at normal incidence (see Fig. [8] for a schematic diagram of the incident polar- 
ization), which have an effective dimension of one. For simplicity, we take the surrounding 
medium to be air. Although, it would be straightforward to introduce other dielectric layers 
into the calculations. Because wavevector components only exist for the direction normal 
to the surface, these systems are ideal for studying and qualitatively highlighting nonlocal 
effects. Furthermore, they allow us to draw some connections with related experimental 
results 20 . 

The transmission, reflection, and absorption spectra for 2, 10, and 20 nm thick films are 
shown in Figs.[5]-[7j respectively. In the absorption spectra, narrow additional (anomalous) 
absorption peaks occur in the nonlocal results, relative to the local ones. The appearance 
of these peaks is identical to theoretical predictions^ and experimental observations^ on 
other thin metal films, where they are the result of optically excited longitudinal (or vol- 
ume) plasmons. (The effects are called such because they are longitudinal to k and are 
contained within the volume of the structure, unlike surface plasmons, which propagate 
along a metal-dielectric interface.) Not surprisingly, at the anomalous absorption energies, 
there is corresponding decrease in the transmission. However, contrary to the initial expec- 
tation of a decrease in reflection, we find that there can be either an increase or a decrease, 
depending on if the corresponding absorption occurs well above (giving an increase) or below 
(giving a decrease) the surface plasmon energy, which is around 2.65 eV for the 10-nm film, 
for example. 

Although a little hard to discern from Figs. [5] — [T], but can be inferred from previous 
results^, the anomalous absorption peaks redshift as the film thickness is increased. This 
causes many more such peaks that were at higher energies to appear in the optical range. 
For example, there are three peaks for the 2-nm film (Fig. [5|, but twelve for the 10-nm one 
(Fig. [6]). In addition, their intensities drastically decrease, where by 20-nm, the nonlocal 
results are almost converged to the local ones; Fig. [7j We will revisit these points below. 

In order to determine whether the anomalous absorption in these results is actually from 
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FIG. 5: (color online) (top) Transmission, (middle) reflection, and (bottom) absorption spectra for 
a 2 nm thick Au film illuminated at normal incidence. Both local (broken red lines) and nonlocal 
(solid blue lines) calculations are shown. 

the excitation of longitudinal plasmons, intensity profiles of |D| 2 at the anomalous absorption 
energies for the 2-nm film [1.14 (not shown in Fig. [5]), 3.36, and 5.54 eV] can be examined; 
Fig. [8] Well-defined standing- wave patterns of |D| 2 are seen inside the films longitudinal to 
k, confirming the assumption of their nature. The wavelengths of these standing waves is 
found to satisfy the condition 

A L = 2d/m , (16) 
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FIG. 6: (color online) (top) Transmission, (middle) reflection, and (bottom) absorption spectra for 
a 10 nm thick Au film illuminated at normal incidence. Both local (broken red lines) and nonlocal 
(solid blue lines) calculations are shown. 



where d is the film thickness and m 



1,3,5,..., which means that odd numbers of half- 



wavelengths fit longitudinally into the structures. The wavelengths defined by Eq. ( 16 ) will 



hereon be referred to as "modes", characterized by m. In Fig. [8j it is the m = 1,3, and 5 
modes that are explicitly shown. These results are significantly different from the local ones, 
where relatively uniform |D| 2 patterns are found at all energies (not shown). Nonetheless, 
as mentioned above, this analysis agrees with previous results on analogous systems^^, 



15 



0.35 

S 

o 

SB 

.3 0.25 



£ 0.15 



0.75 

_o 

"5 0.55 
u 
-S 

0.35 



O 



0.45 



0.35 



O 

< 



£ 0.25 



0.15 
0.05 



■ local 
■nonlocal 




■ local 
■nonlocal 




--local 
— nonlocal 



1.0 2.0 3,0 4.0 

Energy (eV) 



5.0 



6.0 



FIG. 7: (color online) (top) Transmission, (middle) reflection, and (bottom) absorption spectra for 
a 20 nm thick Au film illuminated at normal incidence. Both local (broken red lines) and nonlocal 
(solid blue lines) calculations are shown. 

providing further support for the validity of our method. Based on the observations in Fig. 
[8] and the analysis above, it makes sense that the anomalous absorption features should 
redshift with increasing thickness and that their intensity should decrease with increasing 
m. 



From the discussion above, and that in Subsection |II B[ the approximate anomalous 
absorption energies can be predicted. From Eq. ([9]), it is seen that rapid variations in e(k, u) 



16 




0.0 0.5 1.0 

FIG. 8: (color online) Normalized |D| 2 intensity profiles inside a 2 nm thick Au film at energies of 
(left) 1.14, (middle) 3.36, and (right) 5.54 eV. The polarization and direction of incident light is 
indicated; in each image, the sides of the metal film are padded on the left and right using solid 
black lines. 

will occur when u ~ (3k, which will likely lead to an absorption condition; and from Eq. 



(16), it is seen that longitudinal plasmons with wavelength Al are excited inside a structure 
of thickness d, which will result in momentum states with magnitude k = Thus, 
everything needed to predict the approximate anomalous absorption energies is known: Tluj 
= mfin /d. Using the 2-nm film as an example, this analysis predicts anomalous absorption 
at energies of hu = m ■ 1.44 eV. For the first three m modes, these are 1.44, 4.31, and 7.19 
eV, while those rigorously calculated are 1.14, 3.36, and 5.54 eV, respectively; Fig. |5j While 
not exact, the predictions are reasonably close. Part of these differences can be attributed to 



the grid spacing error, as outlined in Section IV , which in this case leads to an uncertainty in 



d of ±0.2 nm. This analysis could also be applied to related experimental results^. However, 
it is important to keep in mind that this is a simple approximation, and if appropriate, more 
accurate values should be obtained from full calculations or, if possible, by using rigorous 
theoryQS. 

B. Spherical Nanoparticles (3D Systems) 

In this subsection, we study spherical nanoparticles, utilizing the full 3D nonlocal electro- 
dynamics method outlined in Subsection |II A and the Appendix. The optical responses of 



nanoparticles with diameters of 4, 7, and 15 nm are shown in Fig. M (The polarization and 
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FIG. 9: (color online) Extinction cross sections of spherical Au nanoparticles with diameters of 
(top) 4, (middle) 7, and (bottom) 15 nm. Both local (broken red lines) and nonlocal (solid blue 
lines) calculations are shown. 

direction of incident light is irrelevant.) Note that for small nanoparticles, such as these, 
the optical responses are predominately absorption, as scattering does not play a significant 
role for sizes less than approximately 20 nm. Figure [9] shows that inclusion of nonlocal 
effects results in significant anomalous absorption and LSPR blueshifting for both the 4 and 
7-nm nanoparticles. In fact, these effects are so large that the main LSPRs are hardly even 
distinguishable. 
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As discussed in Subsection V A the anomalous absorption peaks arise from the excitation 
of longitudinal plasmons. However, unlike the systems discussed in that subsection, this 
effect diminishes much faster as the nanoparticle size is increased, such that at 15-nm, the 
anomalous peaks show up only as slight indents on the main LSPR. These differences can be 
attributed to two effects. First of all, in a spherical nanoparticle, scattering of the incident 
light off of the exterior surface generates many k-components that can interact and dephase 
one another, especially for the high-order m modes with multiple nodes. Secondly, scattering 
of the conduction electrons that compose the longitudinal plasmons off of the interior surface 
can also lead to dephasing. Each of these processes causes nonlocal effects to diminish at 
much smaller distances in spherical nanoparticles, relative to metal films. 

LSPR blueshifting is most apparent for the 15-nm nanoparticle. This is simply because 
the anomalous absorption is low, allowing this peak to be clearly identified. The local LSPR 
is seen at 2.57 eV, while the nonlocal one is seen at 2.71 eV. This effect can be understood 
by looking at the form of Eq. ([9]). When nonlocal effects are included, the interplay between 
co and k causes all effects {e.g., the absorption condition of Real[£:(k, u)] = —2, in this case} 
to appear at higher energies compared to k = 0. 

Based on the results in Fig. [9j one might wonder why such strong nonlocal effects have 
not been experimentally observed in such systems. Obviously these effects are important, 
and have been observed in other cases^. There are many possible reasons for this. The most 
probable one is that experimental measurements are often made on heterogeneous collections 
of nanoparticles. Given that nonlocal effects are very sensitive to nanoparticle dimensions 



(see Section IV, for example), slight heterogeneity could effectively average them away. 
Support for this claim comes from an experimental study of isolated Au nanoparticles^, which 
clearly demonstrated the LSPR blueshift, and possibly anomalous absorption features^. 
Another possible explanation is that our choice of f3 2 is not optimal (which is directly related 
to the strength of nonlocal effects), as we have recently argued for metallic nanoshells^. The 
hydrodynamic Drude model neglects quantum mechanical exchange and correlation effects, 
which in a local density approximation would decrease 1 . A third possible explanation 
is that our choice of damping parameter A is too low. Increasing this would smooth all 
spectral features (i.e., the anomalous absorption would not appear as strong). Support for 
this comes from a combined theoretical (local electrodynamics) and experimental study of 
metallic nanoshells, where values of greater than 1.0 are needed to describe the results (this 
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corresponds to L e g reduced below that based on geometric considerations alone )PH 
C. Nanowires (2D Systems) 

In our previous Letter^, we demonstrated that nonlocal effects are particularly important 
in structures with apex features, such as triangular nanowires. In such structures, optical 
responses can be affected by nonlocal effects for much larger sizes than for those with smooth 
geometries, such as cylindrical nanowires. Additionally, even though far-field optical prop- 
erties seem to converge to local ones at large sizes with regard to anomalous absorption 
and LSPR blueshifting (to a large extent), the near-field properties, such as electric field 
enhancements, hereon referred to as |E| 2 enhancements, do not. 

As a final application of the method presented in Subsection II A| and the Appendix, we 



study the optical responses and |E| enhancements around isolated cylindrical, square, and 
triangular nanowires with 50 nm diameters or side-lengths, common sizes used in experi- 
mental and theoretical studies. |E| 2 enhancements from such structures have been studied 
in the past! 36 * 37 ! However, to the best of our knowledge, all previous studies have been car- 
ried out using local electrodynamics (at least for non-cylindrical structures), except for the 
aforementioned discussion in our previous Letter^. 

Optical responses of nanowires illuminated with the E-field polarized along the long axis 



of each structure are shown in Fig. 10 (Schematic diagrams of the polarization are shown 
in Figs. 11 and 12 ) Because of their large sizes, these structures do not exhibit distinct 
anomalous absorption. Nonetheless, many closely spaced longitudinal plasmon modes do 
exist (vide infra), which leads to very minor, closely spaced "bumps" in the nonlocal results, 
as well as LSPR blueshifting. These modes can again be confirmed by looking at intensity 
profiles of |D| 2 ; Fig. 11 Unlike the results in Fig.pl the longitudinal plasmons in Fig. 11 form 



much more complex patterns. These can be attributed to two (related) effects. One is that 
the size of the structure along the longitudinal direction of the incident field is not the same 
at all positions. Therefore, for a given energy, modes of different order will be sustained at 
multiple positions along the structure, at each place where Eq. ( 16 ) is satisfied^. This is also 



one of the reasons why nonlocal effects are so strong in structures with apex features, and 
why they can remain important in them for arbitrarily large sizes. In other words, low-order 
longitudinal plasmon modes can always be sustained near the apex. And two, scattering 
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FIG. 10: (color online) Extinction cross sections of (top) cylindrical, (middle) square, and (bottom) 
triangular Au nanowires with diameters or side-lengths of 50 nm. Each was illuminated with the 
E-field polarized along the longest axis of the structure. Both local (broken red lines) and nonlocal 
(solid blue lines) calculations are shown. 

of the incident field off of a curved nanowire surface generates many k-components, which 
can excite longitudinal plasmons along directions other than that of the incident k, creating 
an interference pattern. This effect also leads to the dephasing of longitudinal plasmons, as 
discussed in Subsection IV Bl 

It is also interesting to look at intensity profiles of |E| 2 at the LSPR energies, near 
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FIG. 11: (color online) Normalized |D| 2 intensity profiles at the LSPR energies in and around (left) 
cylindrical, (middle) square, and (right) triangular Au nanowires with diameters or side-lengths 
of 50 nm. Both (top) local and (bottom) nonlocal calculations are shown. The polarization and 
direction of incident light is indicated; the nanowires are outlined in white. 



where this quantity is expected to be maximized^; Fig. 12 (Note that these energies are 
slightly different in the local and nonlocal results, due to LSPR blueshifting.) In Fig. 12, the 
|E| 2 profiles have been normalized for each geometry so that relative comparisons between 
the local and nonlocal results can be made. Because of this, it is not possible to directly 
compare the results for different geometries, but such comparisons will be considered below. 
In all cases, the |E| 2 values are qualitatively similar, both inside and around the structures. 
Quantitatively, however, the nonlocal enhancements are clearly lower. This is especially true 
for the triangular nanowires. 

In order to accurately assess the |E| 2 enhancements around the nanowires, the precise 
maximum and average values were determined; Table [TJ The average values refer to fields 
averaged over certain distances from the nanowire surfaces. In all cases, decreases in both 



quantities are seen in the nonlocal results (as could also be inferred from Fig. 12, and the 
discussion above). However, for the cylindrical nanowires, these decreases are negligible. It 
is also interesting to note that the average enhancements are higher 1.0 nm away from the 
surface than they are at 0.5 nm. For the square nanowires, the decreases are noticeably 
larger. There is approximately a 10% difference in both quantities at 1.0 nm, and a 6% 
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FIG. 12: (color online) Normalized |E| 2 intensity profiles at the LSPR energies in and around (left) 
cylindrical, (middle) square, and (right) triangular Au nanowires with diameters or side-lengths 
of 50 nm. Both (top) local and (bottom) nonlocal calculations are shown. The polarization and 
direction of incident light is indicated; the nanowires are outlined in white. 

difference in the average values at 2.0 nm. The decrease in the difference between average 
enhancements at a further distance from the surface is expected, as the near-fields that 
contribute to this exponentially decay. The decreases in |E| 2 enhancements for the triangular 
nanowires are strikingly larger than for the other geometries. For example, decreases of 104% 
and 61% in the maximum and average values, respectively, are seen at 0.5 nm. 

Considering that some physical processes are dependent on |E| 4 enhancements^, such 
as surface-enhanced Raman scattering (SERS), the differences between local and nonlocal 
electrodynamics could have significant implications for the interpretations of results. This 
statement is based on the fact that the nonlocal calculations are, in principle, more rigorous 
than the local ones. For example, if the actual electromagnetic contribution to SERS is 
smaller than expected on the basis of local theory, it is possible that chemical effects play 
a more important role than has been considered in the pastPP. Such results are also likely 
to play a large role in the accurate interpretation of electron energy loss measurements 
for anisotropic nanoparticle structures, which have recently received attention within the 
framework of local electrodynamics®'. 
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TABLE I: Maximum and average |E| 2 enhancements at the LSPR energies around cylindrical, 
square, and triangular nanowires with diameters or side-lengths of 50 nm. Each was illuminated 
with the E-field polarized along the longest axis of the structure. Distances from the nanowire 
surfaces over which averages were obtained are specified. 



Nanowire Shape 


Max. 


Avg. at 0.5 nm 


Avg. at 1.0 nm 


Avg. at 2.0 nm 


Cylindrical (local) 


8.64 


2.42 


2.47 


2.40 


Cylindrical (nonlocal) 


7.85 


2.32 


2.39 


2.34 


Square (local) 


60.58 


3.54 


3.33 


3.01 


Square (nonlocal) 


39.79 


3.02 


2.91 


2.69 


Triangular (local) 


145.77 


5.49 


4.90 


4.18 


Triangular (nonlocal) 


71.40 


3.42 


3.30 


3.01 
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VI. SUMMARY AND OUTLOOK 



In summary, we detailed our electrodynamics method to calculate the optical response 
of an arbitrarily shaped structure described by a spatially nonlocal dielectric function. The 
formulation was based on converting the hydrodynamic Drude model into an equation of mo- 
tion for the conduction electons, which then served as a current field in the Maxwell- Ampere 
law. By discretizing this equation using standard finite-difference techniques, we incorpo- 
rated it into a self-consistent computational scheme along with the standard equations used 
in the finite-difference time-domain method. 

Using the example of a cylindrical Au nanowire studied in our previous Letter^, we 
demonstrated the remarkable accuracy of our method through comparisons to analytical 
results. As new applications, we calculated the optical responses of thin metal Au films, 
Au nanowires, and spherical Au nanoparticles. These calculations demonstrated a num- 
ber of effects that result from the spatial nonlocality in the dielectric response, including 
anomalous absorption, blueshifting of localized surface plasmon resonances, and decreases 
in electromagnetic field enhancements. 

The results presented demonstrate the importance of including nonlocal effects when 
describing metal-light interactions at the nanometer length scale. It is presently difficult to 
compare our results directly with existing experimental studies, because most of these have 
involved heterogeneous collections of particles or non-continuous systems, which in the small 
size limit tend to average over nonlocal effects. It is our hope that these results will motivate 
new, and more precise experimental studies, particularly those on isolated nanostructures, 
where nonlocal effects are likely to play a large role. 

In the future, we plan to derive more accurate expressions than the hydrodynamic Drude 
model. By incorporating exchange and correlation effects, we also plan to compare nonlocal 
calculations directly with quantum mechanical approaches, such as electronic structure the- 
ory. Such expressions will allow more even more accurate descriptions of nonlocal optical 
phenomena than were presented in this work. 
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Appendix: Nonlocal Finite-Difference Equations 



In this appendix, the finite-difference equations used to model nonlocal dielectric effects 
are derived, and their implementation is discussed. 

First, the temporal derivatives in Eqs. Q and (15) are discretized using a leapfrog 
algorithmpl, 

(A.l) 



Ho^ 2 — — = -V x E(x) n 



At 

^oo E(x) " + ^ E(x)n + J l,(x)" +1/2 + Jhd(x)" +1 / 2 = V x H(xr+V2 (A.2) 



where the superscript n denotes a discrete time step. Equations (13) and (14) are discretized 



using central finite-differences (necessary because of the second-order derivatives) centered 
at time-step n, 



J Lj (x)" +1 - 2J Lj -(x)" + Jl^x)"- 1 
At 2 
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J L ,(x)™ +1 - J Lj (x 



n-l 



2At 



soAe-LjUjlj'- 



^Jl,(x)" = 
E(x) n+1 - E(x 



in— 1 



2At 



(A.3) 



JHD(x)" +1 -2J H D(x)™ + JHD(x) n - 1 , _ J H d(x)" +1 - Jhd(x)^ 1 



At 2 



+ 7" 



2At 



-/3 2 V 2 J HD (x)' 

E(x) n+1 - E(x)™- 1 
2At 



• (A.4) 



Next, update equations for Jlj( x ) and Jhd( x ) are obtained by rearranging Eqs. (A.3) 



and (A.4) 



J L ;(x) n+1 = a Lj J Lj (x) n + ^Jl^x)"- 1 + VLj 



E(x 



,71+1 



E(x 



,71-1 



2At 



(A.5) 



where 



and 



2 - u 2 Lj At 2 
1 + 6 Lj At 
1 - 5 Lj At 
~ ~ 1 + 6 Lj At 

SoAsLjUljAt 2 

1 + 6 Lj At 



Jhd(x 



.71+1 



«HDJHD(x) n + ^HDJHD(x) n + ?7HD 



E(x 



,77+1 



E(x 



,71-1 



2At 



(A.6) 
(A.7) 
(A.8) 

(A.9) 
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where 



4 + 2At 2 /3 2 V 2 
2 + 7At 
2-7At 
~2 + 7 At 

2e WDAt 2 
~ 2 + 7At 



«HD 



(A.10) 
(All) 
(A.12) 



Note that g;hd is an operator, rather than a simple coefficient. To use Eqs. (A. 5) and (A. 9) 



in Eq. (A.2), Jlj( x ) and J H d( x ) ar e centered at time-step n + 1/2 by averaging, 



1/2 _ Jl,(x)" +1 + J^(x)' 



Jnnfx)" +1 / 2 



Jhd(x)™ +1 + Jhd(x)™ 



(A.13) 



(A.14) 



Equations (A.2), (A. 5), and (A.9) all contain E(x) n+1 . To obtain a consistent update, 
Eqs. (A.13) and (A.14) [using Eqs. (A. 5) and ( A.9[ )] are inserted into Eq. (A.2) and rear- 
ranged, 

1 



Efx 



,n+l 



Cl+C2 



CiE(x) n + C 2 E(x) n - 1 + V x H(x) n+1/2 - J T (x) n ' n " 1 



where 



EoE c 



At 



(A.15) 

(A.16) 
(A.17) 



and 



Jt(x) 



n,n— 1 



He 



K i + l)J Li (xr + ^JL,(x)"- 1 

(anD + ^JuD^T + ^DJuB^Y 1 - 1 } ■ (A.18) 



Rearrangement of Eq. (A.l) gives the appropriate update equation for H(x) 



H (x)™ +1 /2 = H(x)™- 1 / 2 - —V x E(x 



(A.19) 



In order to satisfy Eqs. ^ and ([6]), a Yee spatial discretization^ is used for the compo- 
nents of E(x) and H(x) (i.e., they are offset and circulate one another). The Jlj(x) and 
Jhd(x) components are centered at the same spatial locations as the corresponding E(x) 
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components. All of the spatial derivatives in Eqs. (A.9), (A. 15), (A. 18), (A. 19) and are 
approximated using central finite-differences. 

In order to model an arbitrarily shaped structure, the Jlj(x) and Jhd( x ) components 
only exist at the grid positions of the corresponding nonlocal material. By not updating the 
currents outside of the structure, the additional boundary condition of Pekar is imposed^ 
- i.e., the total nonlocal polarization current vanishes outside of the structure. 



Equations (A. 5), (A.9), (A. 15), and (A. 19) form the complete and consistent set necessary 
to solve Eqs. (j3]) - ^ for materials described by the constitutive relationship in Eq. ^ with 
the dielectric function given in Eqs. ([9]) . 
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